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Abstract 

We present a numerically stable Quantum Monte Carlo algorithm to calculate 
zero-temperature imaginary-time Green functions G{f, r) for Hubbard type 
models. We illustrate the efficiency of the algorithm by calculating the on- 
site Green function G{f = 0, r) on 4 x 4 to 12 x 12 lattices for the two- 
dimensional half-filled repulsive Hubbard model at U/t = 4. By fitting the 
tail of G{f = 0,r) at long imaginary time to the form e~'^^'=, we obtain a 
precise estimate of the charge gap: Ac = 0.67 it 0.02 in units of the hopping 
matrix element. We argue that the algorithm provides a powerful tool to 
study the metal-insulator transition from the insulator side. 
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I. INTRODUCTION 



The purpose of this article is to describe a numerically stable Quantum Monte Carlo 
(QMC) algorithm to calculate zero-temperature imaginary time displaced Green functions: 
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where c,(r) = e'^^-'^^^c^e-'^^"^^). (1) 

Here |\l/o) denotes the ground state of the considered Hamiltonian H, c[. creates an electron 
with quantum numbers x, 0(t) is the Heaviside function and n is the chemical potential 
which has to satisfy: 

Tr (e-/^(^-^^)iV) ^ (^oliVl^,) 
i"" Tr(e-/^(^-''^)) (^ol^o) ' ^ ' 

in the zero-temperature limit T = — 0. The above equation implies that the chemical 
potential corresponding to the desired particle density, has to be known prior to the 

simulation. In a metallic state n(/i) is in general not known a priori. However, in an 
insulating state at zero-temperature n(/i) is constant and in general known for chemical 
potentials within the charge gap Ac. In this situation, the here described algorithm proves 
to be a powerful tool. The above T = Green functions have already been calculated with 
QMC methods by Deisz et al Since their algorithm does not incorporate a numerical 
stabilization scheme, they are restricted to relatively small values of r (i.e. r = 2.5 in units of 
the hopping matrix element for the one-dimensional Hubbard model). This articles follows 
the work of Deisz et al. and describes a numerical stabilization scheme which allows one to 
calculate T = Green functions for arbitrary values of r. 

To demonstrate the efficiency of the algorithm, we calculate Gx,y{T) for the two- 
dimensional half-filled Hubbard model: 

^ = E iT^M. + (n,, - I) (n. , - 1) . (3) 
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The quantum numbers i and a denote lattice site and ^-component of spin respectively, 
= ct^c^^ , and T^j = —t li i and j are nearest-neighbors. In this notation half-band 
filling corresponds to /i = 0. As a non-trivial test of the algorithm, one may fit the tail of 
Gx,x{j) = {i,cr)) to the form e~'^^'= to obtain the charge gap Ac. At U/t = 4 and after 
extrapolation to the thermodynamic limit, our QMC data yields Ac/t = 0.67 ± 0.02. This 
value stands in good agreement with previously determined values of Ac 0. 

The article is organized in the following way. In the next section, we briefly describe 
the zero-temperature auxiliary- field QMC algorithm for the Hubbard model . We then 
present our solution for the numerical stabilization of the time displaced Green functions. 
In section 3 we describe our calculation of the charge gap for the two-dimensional Hubbard 
model at U/t = 4. In the last section, we draw some conclusions and discuss the potential 
applications of the algorithm. 

II. THE ZERO-TEMPERATURE QMC ALGORITHM 

Since the Hubbard model conserves particle number, and we are working in the canonical 
ensemble, one may factorize the chemical potential to obtain: 



G.,y{T) = e(r) 
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Due to the above relation we consider the calculation of the T = Green functions at = 0. 

The idea behind the zero temperature QMC algorithm is to filter out the ground state 
from a trial wave function |\I/t) which is required to be non-orthogonal to the ground state: 
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The QMC calculation of 

proceeds in the following way. The first step is to carry out a Trotter decomposition of the 
imaginary time propagation: 
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Here, Ht (Hu) denotes the kinetic (potential) term of the Hubbard model and mAr = 2G. 
Having isolated the two-body interaction term, Hu, one may carry out a discrete Hubbard 
Stratonovitch (HS) transformation |^ to obtain: 

^-ArHu = c E exp f E ciD^A^)Cy] , (8) 
s \x,y / 

where s denotes a vector of HS Ising fields. For the Hubbard model (|^), we take: 

Dr.jA^ = S^yS.jcosh~\ATU/2)s.a. (9) 

The constant C = exp(— ArA^f//2)/2^ for the A^-site system will be dropped below. The 
imaginary time propagation may now be written as: 

e-'®^ = E^K2e,0) + O((Ar)2) 

s 

where C/,-(2G, 0) = e'^^^^/^ fj^ 4d,,,(4)c,) e-^^^*/l (^q) 

n=l \x,y J 

In the above notation, G^y{Q, r) is given by: 

The trial wave function is required to be a Slater determinant: 

I^T)=nfE4^x,n)|0). (12) 

n=l \ X J 

Here A''^ denotes the number of particles and P is an Ng x Np rectangular matrix where A^^ 
is the number of single particle states. Since Ug{2Q, 0) describes the propagation of non- 
interacting electrons in an external HS field, one may integrate out the fermionic degrees of 
freedom to obtain: 

^^.(e, r) = X: P. + r, 6)) G,(e, 6)],.,,^ + 0((Ar)2). (13) 

s 

In the above equation. 
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Bg{Q2, 0i) = n e-^^^/2e-°(^"^e-^^^/2 ^^ere mAr = Oi and riaAr = 62, 

n=ni+l 



Mj=P^S,-(2e,0)P, 
det(Mj) 
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and (G,(e,e)),,^ = (^,1^,(29,0)1^.) = 

(7 - B,(e, o)PMi^p^B^(2e, e)) . 

Here I is the unit matrix, I^^y — 5x,y In the same notation one obtains: 

= Y.Ps\{Gs{Q,Q)-I)Bz\Q + tM , t>0. (14) 

^ L lx,y 
s 

Summarizing, the zero-temperature imaginary-time Green function may be calculated from: 

Gx,y{T) = Jim (0(r)G>,(0, r) + e(-r)G<,(e, -r)) + 0((Ar)^). (15) 

At half-band filling and due to particle hole symmetry, one may chose a trial wave 
function such that Pg is positive definite. Pg may be interpreted as a probability distribution 
and sampled with Monte-Carlo methods. 

A. Numerical Stabilization 

The origin of the numerical instabilities occurring in the calculation of Green functions 
may be understood by considering free electrons on a two-dimensional square lattice. 



<?J> 

Here, the sum runs over nearest-neighbors. For this Hamiltonian one has: 

(*ol4WCikl*o) - exp (T(e,- - ^J)) (*ol4cfel*o), (17) 



where eg = —2t{cos{kax) + cos{kay)), a^, Sy being the lattice constants. The chemical 

potential satisfies equation (|^) and we will assume |\E'o) to be non-degenerate. In a numerical 

calculation the eigenvalues and eigenvectors of the above Hamiltonian will be known up to 

machine precision, e. In the case eg — > 0, (\&o|cgCg|^'o) = 0. However, on a finite 

precision machine the later quantity will take a value of the order of e. When calculating 

(\E'o|ct(r)cg|\E'o) this roundoff error will be blown up exponentially and the result for large 

values of r will be unreliable. In order to circumvent this problem, one may do the calculation 

at finite temperature and then take the limit of vanishingly small temperatures: 

, exp ( rfer — u)] 
(^ol4(r)cg|vI/o) = hm ^ (18) 



Even if the eigenvalues are known only up to machine precision, the right hand side of the 
above equation for large but finite values of /5 is a numerically stable operation. Although 
very simple, this example refiects the underlying numerical instabilities occurring in the 
calculation of the Green functions. 
We now consider the calculation of 

Gg{e + r, 0) = Bg{e + r, 0)G,-(e, 6) and 

G,{Q, e + r) = (a-(e, e) - /) B^\e + r, e) (19) 

required to compute G^y{Q, r) (see equation (|13D) and ^^^^(6, r) (see equation (|I4D) respec- 
tively. The equal-time Green functions, Gg{Q,Q), may be calculated to machine precision 
1^-^. The matrices Bg{Q + T, 0) contain scales which grow and decrease exponentially with 
r. As in the above example, a straightforward multiplication of both matrices will lead to 
numerical instabilities for large values of r. Here, the problem is much more severe since 
the presence of the HS field mixes different scales. In order to circumvent this problem, we 
propose the following stabilization scheme. 

Since the trial wave function is a Slater determinant, we can find a single particle Hamil- 
tonian, Hq = J2x,y(aho)x^yCy, which has |\I't') as a non-degenerate ground state. The equal 
time Green functions may then be written as: 
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Tr(e-/5^^of/.(2e,0)) /3-oo 

The last equality follows after integration of the fermionic degrees of freedom. Inspiring 
ourselves from the work of Hirsch P] we calculate the time displaced Green functions in 
equation (p!9|) with: 
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For very large but finite values of /3, we can calculate the left hand side of the above equation 
by using matrix stabilization techniques developed for finite temperature QMC algorithms. 
The basic idea behind those numerical stabilization techniques is to keep the different scales 
occurring in the matrices Bg separate (for a review see reference [^J). This is achieved by 
decomposing the matrices Bg into a UDV form where U is an orthogonal matrix, D a 
diagonal matrix containing the exponentially large and exponentially small scales, and V a 
triangular matrix. The calculation of the left hand side of the above equation is done in the 
following way: 
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In the above equation, the matrix D3 contains only exponentially large scales since the 
matrices (V2?7i)^^ and act as a cutoff to the exponentially small scales in the 

matrices D2 and Di. Since the other matrices are all well conditioned, the final matrix 
multiplication is well defined. 

A convenient choice of Hq is obtained in a basis where the trial wave function may be 
written as: 

I*t) = 117110). (23) 

n=l 

In this basis, we define Hq through 

^,-,110) = ( -'il°> = (24) 

[ +7*10) if 7hv.|*r) =0 

(Here, the energy unit is set by the hopping matrix element t.) For this choice of values 
of pt ~ 40 were well sufficient to satisfy equation (^O]) within required numerical precision 
||1CI|| . The above numerical stabilization scheme was indeed successful in all examined cases. 

III. EVALUATION OF THE CHARGE GAP FOR THE TWO DIMENSIONAL 

HUBBARD MODEL 

We carried out our simulations on 4 x 4 to 12 x 12 lattices for the two-dimensional half- 
filled (/i = 0) repulsive Hubbard model at U/t = 4. Periodic boundary conditions were 
assumed. A spin singlet ground state of the kinetic energy in the Hubbard Hamiltonian was 
used as a trial wave function. We test the quality of this trial wave function on a 6 x 6 
lattice. Figure 1 plots (^r|e~®^Oe^®^|\E'r)/(^I^T|e~^®^|\l'T) as a function of the projection 
parameter 6 for O = 5'(7r,7r)/A^ = g^Z^^exp (iQrj S{f) ■ S{0) (sohd circles in Figure la) 
and O = E/N = H/N - f//4 (sohd circles in Figure lb). Here Q = (7r,7r)/a, S{r) is the 
spin operator on site r and N denotes the number of sites. Already for values of Qt = 2.5, 
both considered observables have converged within our estimated statistical uncertainty. For 
comparison, we have plotted Tr (^e^^^^o) /Tr (e^^®-^) for the same observables (triangles 
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in Figure 1). Values of Qt at least twice as large are required to obtain approximate ground 
state results. 

Another source of systematic errors comes from the discretization of the imaginary time 
propagation. In Table | the At dependence of the energy and ^(Tr, vr) is given. The data 
are obtained form the zero-temperature QMC algorithm on a 6 x 6 lattice and at 2Qt = 5. 
The values at At = are obtained from a least square fit of the finite At results to the 
form a + 6(Ar)^. We carried out our simulations at Art = 0.125. As may be seen from 
Table |, this value of Art produces a systematic error contained in the quoted errorbars for 
the energy and a systematic error of less than 1% for ^(Tr, vr). 

To obtain an estimate of the charge gap, we consider 

G(r = 0,r) = l^G.,.(r), r > 0. (25) 

Here, x stands for spin and site indices. Inserting a complete set of eigenstates of the 
Hamiltonian H in the + 1 particle Hilbert space yields: 

Gir = 0,r) = l(^o'^|cx|^r^)rexp (-r (e^^ - E^^)) . (26) 

^ n,x 

where H\^^+^) = E^+^\^^+^) and H\^l^) = ^^|^^). At large values of rt, G{f = 0, r) - 
exp (— rAc) where Ac = Eq'^'^ — Eq corresponds to the charge gap. 

Our results are plotted in Figure 2. For those simulations we have chosen Qt = 13.5. 
Since values of r up to T^axt = 12 were considered, the effective projection parameter is 
given by: Qeff = © — Tmaxf^ = 7.5/t. As may be seen from Figure 1, this value of the 
projection parameter is more than sufficient to filter out the ground state from the trial 
wave function. The solid lines in Figure 2 correspond to least square fits of the tail of 



G{f = 0,r) to the form exp(— rA^) [Tl[. The estimated value of the gap as a function of 
system size is plotted in Figure 3. A least square fit of the data to the form a + b/L, where 
L denotes the linear length of the lattice, yields A^/t = 0.67 ± 0.02 in the thermodynamic 
limit. This value stands in good agreement with the value of the charge gap obtained by 
Furukawa and Imada 0: Ac/t = 0.58 ± 0.08 (solid circle in Figure 3.) As may be seen from 



the comparison of errorbars, the accuracy of the estimation has been much improved in the 
present study. 



IV. CONCLUSIONS 

We have presented an efficient, numerically stable, QMC algorithm to calculate T = 
imaginary time Green functions for Hubbard type models. As a non-trivial test application 
of this algorithm, we have obtained an accurate estimate of the charge gap for the two- 
dimensional half-filled repulsive Hubbard model at U /t = 4: A^/t = 0.67 ± 0.02. 

The algorithm is formulated in the canonical ensemble. Hence, the relation n{^) has to 
be known prior to the simulation. This renders the algorithm hard to use in a metallic state. 
However, in an insulating state at zero temperature n(/i) is constant, and generally known, 
for chemical potentials within the charge gap. In this situation the here presented algorithm 
proves to be a powerful tool. We illustrate this by considering the two-dimensional Hubbard 
model. Due to Equation (^), it suffices to know the Green functions at = (half- filling) 
to be able to determine them trivially for any other chemical potential within the charge 
gap. At /X = 0, we are not confronted with a sign problem due to particle-hole symmetry 
and the statistical fiuctuations do not blow-up exponentially with growing lattice sizes and 
projection parameters 0. It is however clear from equation (H) that statistical fiuctuations 
will increase (decrease) exponentially with growing positive values of r for yU > (/i < 0). 
In comparison, finite temperature algorithms in the grand-canonical ensemble, are faced 
with a sign problem away from /x = 0. Hence, statistical fiuctuations grow exponentially 
with growing lattice size and inverse temperature. Away from /i = 0, it is thus extremely 
hard to extrapolate any zero temperature result from the finite temperature grand-canonical 
algorithms for large lattice sizes. This renders the here presented algorithm a powerful tool 



for the study of the metal- insulator transition from the insulator side |13 



10 



ACKNO WLED CEMENTS 



F.F. Assaad thanks the JSPS for financial support. The numerical calculations were 
carried out on the Fujitsu VPP500 of the Supercomputer Center of the Institute for Solid 
State Physics, Univ. of Tokyo. This work is supported by a Grant-in-Aid for Scientific 
Research on the Priority Area "Anomalous Metallic State near the Mott Transition" from 
the Ministry of Education, Science and Culture, Japan. 



11 



REFERENCES 

[1] J.J. Deisz, W. von der Linden, R. Preuss and W. Hanke, to appear in Computer sim- 
ulations in Condensed Matter Physics VIII, Eds. D.P. Landau, K.K. Mon and H.B. 
Schiittler (Spinger Verlag, Heidelberg, Berlin, 1995) 

[2] N. Furukawa and M. Imada, J. Phys. Soc. Jpn. 61 (1992) 3331. 

[3] G. Sugiyama and S.E. Koonin, Anals of Phys. 168 (1986) 1. 

[4] S. Sorella, S. Baroni, R. Car, And M. Parrinello, Europhys. Lett. 8 (1989) 663. S. Sorella, 
E. Tosatti, S. Baroni, R. Car, and M. Parinello, Int. J. Mod. Phys. Bl (1989) 993. 

[5] M. Imada and Y. Hatsugai, J. Phys. Soc. Jpn. 58 (1989) 3752. 

[6] P.P. Assaad, Helvetica Pysica Acta 63 (1990) 580. 

[7] J.E.Hirsch, Phys. Rev. B. 28 (1983) 4059. 

[8] J.E. Hirsch, Phys. Rev. B 38 (1988) 12023. 

[9] E. Loh and J. Gubernatis, in Modern Problems of Condensed Matter Physics, edited 
by W. Hanke and Y.V. Kopaev, (North Holland, Amsterdam, 1992), Vol 32, p. 177. 

[10] On average equation ( pOD was satisfied up to 10^^^ in absolute values. 

[11] Since the statistical fiuctuations for different values of r are not independent, we have 
done our error analysis in a basis where the covariance matrix is diagonal (see reference 
0) 

[12] M. Jarrel and J.E. Gubernatis, Bayesian Inference and the Analytic Continuation of 
Imaginary- Time Quantum Monte Carlo Data. Preprint. 

[13] P.P. Assaad and M. Imada, unpublished. 

Figure captions 



12 



Fig. 1 •: (\1/T|e ^^Oe ^^I^^t) / {"^tI^ as a function of the projection param- 

eter O. The trial wave function is a singlet ground state of the kinetic energy in 
Hubbard Hamihonian. A: Tr (g-^e^^o) /Tr (e-^©^). 
We have considered two observables: a)0 = ^(Tr, 7r)/A^, b) O = E/N. 

Fig. 2 lnG(r = 0, r) as a function of lattice size. The solid lines are least square fits 
of the tail of G{f = 0, r) to the form e''^^". 

Fig. 3 Ac as a function of 1/L where L corresponds to the linear size of the square lat- 
tice. The solid circle at 1/L = corresponds to Furukawa and Imada's estimate 
of Ac (see reference 
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TABLES 

TABLE I. Energy per site and ^(Tr, tt) as a function of Art for the Hubbard model aX U/t = A 

on a 6 X 6 lattice. The simulations were carried out with the zero-temperature QMC algorithm at 
2Qt = 5. The quoted values at Art = are obtained from a least square fit to the form a + b{ATt)^ 



Art E/Nt S{7r,7r)/N 

0.0 -0.8575 ± 0.0003 0.1579 ± 0.0007 

0.0625 -0.8571 ± 0.0004 0.1578 ± 0.0009 

0.1 -0.8571 ± 0.0003 0.1570 ± 0.0006 

0.125 -0.8570 ± 0.0003 0.1564 ± 0.0007 

0.166 -0.8563 ± 0.0003 0.1557 ± 0.0008 
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Figure 3 



